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Version  History 


This  release  of  WAVEMAKER  software  incorporates  a  recent  technical  development  achieved 
at  the  Reliability  of  Marine  Structures  Program:  the  ability  to  successfully  predict 
wave  histories  at  other  locations  based  on  the  underlying  first-order  wave  components 
for  a  given  target  observed  wave  history.  WAVEMAKER  3.0  is  not  fully  backward  com¬ 
patible,  that  is,  results  from  a  run  of  earlier  versions  of  WAVEMAKER  can  not  be  exactly 
obtained  with  this  new  release  because  there  have  been  some  minor  improvements 
to  existing  parts  of  the  program.  The  input  files  to  earlier  versions  can,  however,  be 
directly  used  in  this  new  version.  A  version  history  of  WAVEMAKER  follows: 

•  Version  1.0:  Released  in  April  1995  and  documented  in  Jha  and  Winter- 
stein,  1995,  Report  RMS-17,  contains  simulation  capabilities  for  second- 
order  random  waves. 

•  Version  1.1:  Released  in  August  1995,  includes  modification  of  tempo¬ 
rary  intermediate  output  file  to  use  less  disk  space  and  to  reduce  program 
execution  time  by  approximately  50%.  Additionally,  a  DOS  executable  of 
WAVEMAKER  was  included  in  this  release. 

•  Version  2.0:  Released  June  1996,  includes  identification  capabilities  so 
that  underlying  first-order  wave  history  can  be  retrieved  from  an  observed 
wave  history.  Report  RMS-22 

•  Version  3.0;  Released  June  1998,  includes  prediction  capabilities  so  that 
wave  histories  at  other  spatial  locations  can  be  predicted  from  a  user- 
specified  wave  history  at  a  known  location.  Also,  endpoint  discontinuities 
in  the  wave  records  are  treated  prior  to  application  of  FFT’s. 
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Abstract 


WAVEMAKER  is  a  FORTRAN  subroutine  to  simulate  random  non-Gaussian 
ocean  wave  histories.  It  generates  a  first-order  (Gaussian)  wave  process 
with  an  arbitrary  power  spectrum,  and  applies  nonlinear  corrections  based 
on  second-order  hydrodynamics.  Inputs  to  the  routine  include  the  first- 
order  spectrum,  the  water  depth,  and  a  set  of  locations  in  the  along-wave 
direction  at  which  wave  elevation  histories  are  desired.  It  may  thus  provide 
useful  input  to  estimate  loads  on  spatially  distributed  ocean  structures  and 
ships. 

The  WAVEMAKER  package  also  includes  a  separate  driver  program,  which 
facilitates  input/output  and  generates  several  analytical  spectral  models. 
Its  input  is  specified  in  command-line  format,  similar  to  that  of  the  TF- 
POP  program  for  hydrodynamic  post-processing  also  developed  in  the 
Stanford  RMS  program.  An  example  problem  is  included  to  demonstrate 
the  use  of  WAVEMAKER  and  its  driver. 

In  terms  of  methodology,  WAVEMAKER  first  uses  standard  frequency  do¬ 
main  methods  to  generate  first-order  Gaussian  histories  at  each  location. 
For  each  of  these,  WAVEMAKER  then  evaluates  the  full  set  of  second-order 
corrections  according  to  hydrodynamic  theory.  Thus  the  first-order  wave 
process,  with  N  components  at  frequencies  Uny  gives  rise  to  a  total  of  N 
corrections,  spread  over  all  sum  frequencies  a;„  -f  Um,  and  to  another 
corrections  over  all  difference  frequencies  Un  ~ 

WAVEMAKER  also  includes  the  ability  to  identify  the  underlying  first- 
order  Gaussian  history  from  a  given  observed  time  history.  This  feature 
is  particularly  attractive  for  use  in  situations  where  the  second-order  non¬ 
linearity  in  the  waves  is  built-in  into  the  structural  response  calculations. 
To  avoid  double-counting  therefore,  the  input  waves  should  be  filtered  to 
remove  any  second-order  nonlinearity.  WAVEMAKER  takes  in  an  input  wave 
history  and  identifies  its  first-  and  second-order  wave  components.  This 
identification,  an  inverse  feature  to  simulation,  is  based  on  a  Newton- 
Raphson  scheme  to  solve  N  simultaneous  nonlinear  equations  to  identify 
the  first-order  waves  which,  when  run  through  the  second-order  wave  pre¬ 
dictor,  matches  the  observed  waves. 

Wavemaker’s  final  capability  is  to  use  the  identified  underlying  first- 
order  Gaussian  wave  history  to  predict  consistent  first-  and  second-order 
wave  histories  at  alternative  user-specified  spacial  locations. 
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Introduction  to  WAVEMAKER  3.0 


This  release  of  WAVEMAKER  software  incorporates  a  recent  technical  development  achieved 
at  the  Reliability  of  Marine  Structures  Program:  the  ability  to  successfully  predict 
wave  histories  at  other  locations  based  on  the  underlying  first-order  wave  components 
for  a  given  target  observed  wave  history.  WAVEMAKER  3.0  is  not  fully  backward  com¬ 
patible,  that  is,  results  from  a  run  of  earlier  versions  of  WAVEMAKER  can  not  be  exactly 
obtained  with  this  new  release  because  there  have  been  some  minor  improvements 
to  existing  parts  of  the  program.  The  input  files  to  earlier  versions  can,  however,  be 
directly  used  in  this  new  version.  A  version  history  of  WAVEMAKER  follows: 

•  Version  1.0:  Released  in  April  1995  and  documented  in  Jha  and  Winter- 
stein,  1995,  Report  RMS-17,  contains  simulation  capabilities  for  second- 
order  random  waves. 

•  Version  1.1:  Released  in  August  1995,  includes  modification  of  tempo¬ 
rary  intermediate  output  file  to  use  less  disk  space  and  to  reduce  program 
execution  time  by  approximately  50%.  Additionally,  a  DOS  executable  of 
WAVEMAKER  was  included  in  this  release. 

•  Version  2.0:  Released  June  1996,  includes  identification  capabilities  so 
that  underlying  first-order  wave  history  can  be  retrieved  from  an  observed 
wave  history.  Report  RMS-22 

•  Version  3.0:  Released  June  1998,  includes  prediction  capabilities  so  that 
wave  histories  at  other  spatial  locations  can  be  predicted  from  a  user- 
specified  wave  history  at  a  known  location.  Also,  endpoint  discontinuities 
in  the  wave  records  are  treated  prior  to  application  of  FFT’s. 

In  this  report.  Chapter  2  includes  documentation  of  the  simulation  capabilities 
and  is  largely  taken  from  Report  RMS-17.  Chapter  3  documents  the  identification 
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analysis  capabilities,  and  is  largely  taken  from  Report  RMS-22.  Chapter  4  docu¬ 
ments  the  newly  developed  prediction  capabilities  of  Wavemaker.  The  appendices 
include  sample  input  and  output  files  for  the  simulation,  identification  and  prediction 
examples. 


Chapter  2 


Simulation  of  Second- Order  Random 
Waves 


2.1  Introduction 


It  is  common  in  many  ocean  engineering  problems  to  seek  to  simulate  a  time  trace 
of  the  wave  elevation  process,  r){t),  at  one  or  more  locations  in  the  along-wave  direc¬ 
tion.  It  is  most  typical  to  use  a  Gaussian  model  of  77(f)  for  this  simulation,  which  is 
consistent  with  linear  wave  theory.  This  is  due  primarily  to  the  ease  of  simulating 
Gaussian  processes,  e.g.  with  FFT  (Fast  Fourier  Transform)  methods  for  an  arbitrary 
wave  spectrum  (e.g.,  Borgman,  1969). 

We  seek  here  to  demonstrate  and  facilitate  a  similar  frequency-domain  simulation 
capability  for  nonlinear  random  waves  at  a  set  of  spatial  locations  (e.g.,  Figure  2.1). 
These  simulations  split  the  wave  elevation  into  a  random  first-order  (linear)  wave  his¬ 
tory,  771(f),  and  a  corresponding  nonlinear  history  772(f)  which  includes  second-order 
corrections.  FFT  techniques  are  used  to  generate  771(f)  with  an  arbitrary  (first-order) 
wave  spectrum,  Sr,i{uj).  Physical  principles  are  used  to  generate  772(f)  from  771(f), 
based  on  second-order  perturbation  analysis  of  the  underlying  nonlinear  hydrody¬ 
namic  problem.  Thus  if  the  first-order  wave  process  has  TV  components,  at  frequen¬ 
cies  cJn,  V2{t)  includes  second-order  corrections,  spread  over  all  sum  frequencies 
Wn  +  ^m,  and  another  corrections  over  all  difference  frequencies  ujn  —  Wm- 

Note  that  these  second-order  wave  models  are  not  novel;  they  date  back  at  least 
to  the  early  1960s  (e.g.,  Longuet-Higgins,  1963).  More  novel,  however,  is  their  recent 
confirmation  with  respect  to  various  statistics  of  field  measurements  (Marthinsen 
and  Winterstein,  1992;  Vinje  and  Haver,  1994)  and  wave  tank  studies  of  still  more 
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Figure  2.1:  Simulated  wave  time  histories  at  specified  spatial  locations 


severe  seas  (Winterstein  and  Jha,  1995).  Such  second-order  wave  models  also  form 
the  basis  of  state-of-the-art  nonlinear  diffraction  analysis  of  floating  structures  (e.g., 
SWIM,  1995;  WAMIT,  1995).  Note  also  that  the  model  of  772 (t)  used  here  varies 
explicitly  with  water  depth,  as  predicted  by  second-order  theory,  to  reflect  increasing 
nonlinearity  as  we  proceed  to  shallower  water  depths. 


2.2  Methodology 


2.2.1  Underlying  Theory  and  Assumptions 

We  first  consider  771  (t),  the  first-order  wave  elevation,  at  a  specific  reference  location 
(say  x=0).  For  either  frequency-domain  analysis  or  time-domain  simulation,  it  is 
convenient  to  write  771  (t)  as  a  discrete  Fourier  sum  over  positive  frequencies  a;*: 

=  +  =  (2.1) 

*=1  *=1 

To  randomize  Eq.  2.1,  the  phases  Ok  are  taken  to  be  uniformly  distributed,  mu¬ 
tually  independent  of  each  other  and  of  the  amplitudes  Ak-  Furthermore,  we  assign 
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random  amplitudes  with  Rayleigh  distributions,  and  mean-square  value 

E[Al]  =  2S^{uk)du}  =  (Tfc ;  du}  =  uJk-  uJk-i  (2.2) 

Finally,  for  purposes  of  simulation  the  lowest  frequency  interval  dw  is  governed  by 
the  total  period  T  of  the  simulation: 

do;  =  —  (2.3) 

Together,  Eqs.  2.1-2.2  ensure  that  each  of  the  N  frequency  components  in  Eq.  2.1 
is  itself  Gaussian.  We  also  caution  against  the  common  use  of  choosing  deterministic 
amplitudes,  Ak=<7k,  particularly  when  interest  lies  in  preserving  higher  moments  of 
rii{t) — or,  similarly,  the  rms  of  second-order  waves,  loads,  and  responses.  Use  of 
deterministic  amplitudes  can  give  unconservative  estimates;  e.g.,  second-order  rms 
values  that  are  on  average  too  small  (Ude,  1994). 

The  resulting  second-order  wave  at  this  elevation,  772  (t),  is  calculated  from  771  (t) 
as 

772(t)  =  77i(t)  4- A772(t)  (2.4) 

in  which  Ari2{t)  includes  second-order  corrections  at  sums  and  differences  of  all  wave 
frequencies: 


m=l  n=l 


In  general,  the  functions  and  are  known  as  quadratic  transfer  functions 
(QTFs),  evaluated  at  the  frequency  pair  Similar  expressions  arise  in  describ¬ 

ing  loads  and  responses  of  floating  structures;  in  this  case  fff  and  are  calculated 
numerically  from  nonlinear  diffraction  analysis  (e.g.,  WAMIT,  1995).  The  leading 
factor  q  is  included  in  Eq.  2.5  to  alert  readers  to  different  QTF  definitions  in  the 
literature:  various  diffraction  analyses  use  q=l  (WAMIT,  1995)  or  9=1/2  (Molin  and 
Chen,  1990). 

In  predicting  motions  of  floating  structures,  in  view  of  the  relevant  natural  periods 
interest  commonly  lies  with  either  (slow-drift)  or  (springing)  but  not  both.  In 
contrast,  in  the  nonlinear  wave  problem  both  sum  and  diff'erence  frequency  effects  play 
a  potentially  significant  role.  Fortunately,  unlike  QTF  values  found  numerically  from 
numerical  diffraction,  closed-form  expressions  are  available  for  both  the  sum-  and 
difference-frequency  QTFs  for  second-order  waves  (e.g.,  Langley,  1987;  Marthinsen 
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and  Winterstein,  1992).  Including  the  effect  of  a  finite  water  depth  d,  for  example, 
the  sum-frequency  QTF  can  be  written  as 


Kn 


tikmkn  _  ±(,  fi  4.  ,,,2  _L  ,,,  \  4-  g 

UmUtn  2g  ^  ^m^n)  +  2 


1 


tanh(A,n  +  K)d 


(2.6) 


in  which  the  wave  numbers  are  related  to  the  frequencies  a;„  by  the  linear  dispersion 
relation.  Note  that  this  QTF  definition  assumes  g=l/2  in  Eq.  2.5.  This  is  the 
convention  assumed  in  WAVEMAKER.  The  corresponding  difference-frequency  transfer 
function,  is  found  by  replacing  Un  by  -a;„  in  Eq.  2.6. 


2.2.2  Implementation 

On  input  the  simulation  method  requests  the  desired  number  of  simulated  points, 
npts,  and  the  total  duration  T  to  be  simulated.  To  take  advantage  of  discrete  FFT 
(Fast  Fourier  Transform)  techniques,  it  assumes  a  regular  spacing  dt=T/npts  between 
points.  Eq.  2.1  is  then  rewritten  as 

nptsl2  npts 

ni{t)  =  E  cos(a;*i  -I-  Bk)  =  Re  X^e'*"**  (2.7) 

fc=l  k=\ 

Here  the  Xk  are  complex  Fourier  coefficients.  The  lower  half  of  these  directly 
reflect  both  the  random  amplitude  Ak  and  phase  9k  at  frequency  Uk^k  •  dw: 

Xk  =  ;  k  =  l...npts/2  (2.8) 

The  upper  half  are  in  turn  taken  as  the  complex  conjugates  (the  symbol  “*”)  of  the 
lower  half: 

Xnpts-k  =  X*,;  k==l...npts/2  (2.9) 

This  reflects  that  unique  information  is  contained  only  the  lower-half  frequencies; 
indeed,  any  information  in  the  upper  half  frequencies  (above  the  Nyquist)  is  obscured 
by  aliasing. 

Thus,  the  first-order  wave  process  is  generated  by  assigning  random  Ak  and  6k, 
defining  Xk  from  Eqs.  2.8-2.9,  and  finally  taking  the  inverse  Fourier  transform  to 
recover  the  discretized  time  history  To  see  this,  note  that  since  dt'diJ=2ir /npts, 
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Eq.  2.7  can  be  evaluated  at  t=tj  to  give 

npts 

Viitj)  =  Re  E  (2.10) 

k=l 

This  is  precisely  the  definition  of  the  discrete  FFT. 

As  a  minor  technical  issue,  note  that  the  conjugate  symmetry  here  ensures  that 
the  “Real  Part”  operation  in  Eq.  2.1  is  superfluous;  i.e.,  no  imaginary  component  is 
generated.  Also,  to  conform  with  the  FFT  routine  used  the  array  of  Xk  values  is 
shifted  by  one  index:  i.e.,  Xi  corresponds  to  the  frequency  zero  (the  steady  term, 
defined  as  zero),  X2  to  frequency  Ui=dw,  Xz  to  frequency  0)2=2  •  du),  and  so  forth. 

The  second-order  correction  is  generated  similarly.  Starting  with  Eq.  2.5,  substi¬ 
tuting  q=lf2,  N=npts/2  (Eq.  2.7)  and  Xk  from  Eq.  2.8; 

npts  (2  npts  f  2 

A%(t)=2R*l;  +  (2.11) 

m=l  n=l 

The  leading  factor  reflects  the  product  of  g=l/2  and  a  net  factor  of  4  (since  AnA^  is 
4  lATojATnl).  The  program  then  seeks  to  rewrite  both  the  sum  and  difference  frequency 
contributions  in  a  Fourier  sum  analogous  to  Eq.  2.10.  For  example,  the  sum-frequency 
is  assumed  of  the  form 

npts 

A77+  {tj)  =  Re  (2.12) 

fc=i 

The  output  Fourier  coefficients,  Yk,  are  evaluated  by  equating  Eqs.  2.11  and  2.12.  This 
implies  a  sum  over  all  wave  frequency  pairs  {wmi  (^n)  in  Eq.  2.11  that  give  rise  to  output 
sum  frequency  o)k.  The  difference  frequency  Fourier  coefficients  are  constructed  in  a 
similar  way,  and  added  on  the  sum  frequency  Yk  coefficients.  Once  these  combined 
1*  coefficients  are  found  a  (one-dimensional)  inverse  FFT  is  performed  to  recover  the 
second-order  time  history. 


2.2.3  Multiple  Spatial  Locations 


The  linear  dispersion  relation  can  be  used  to  generalize  Eq.  2.4,  which  generates  a 
first-order  wave  at  reference  location  x—0,  to  any  other  spatial  location  x  in  the  along- 
wave  direction.  The  linear  dispersion  relation  is  first  used  to  find  the  wave  number 
kn  associated  with  each  frequency  a;„  in  Eq.  2.4.  The  modified  first-order  simulation 
then  merely  replaces  Unt  +  On  in  Eq.  2.4  by  -  knX  +  0„.  Equivalently,  the  original 
phases  are  first  shifted  to  before  applying  Eq.  2.8  to  define  the  wave 

Fourier  amplitude  Xn-  These  appropriately  modified  Xn  are  also  used  in  Eq.  2.11  to 
find  the  corresponding  second-order  correction  at  this  new  location. 
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2.3  Input  Specification 


This  section  describes  the  various  inputs  required  by  the  program  and  the  syntax  of 
the  input  to  the  driver  routine  for  WAVEMAKER.  This  input  is  provided  in  the  following 
format: 

keyword  args 

where  keyword  is  a  reserved  word  and  args  are  its  arguments.  A  typical  input  file  is 
in  the  following  format: 


#  Typical  input  file:  syntax  description 

simulate  duration  npts  seed 
depth  value 

psd  psdtype  psd-parameters 

define  varlimit  value 

define  gravity  value 

define  omgmax  value 

write  history  fUenamel  filenames 

write  statistics  filenames  filename4 

location  nloc 

valuel 

values 

valuenloc 


Each  of  these  lines  in  a  typical  input  file  is  explained  below: 

#  Typical  input  file:  S3Titax  description 

Any  line  beginning  with  a  “#”  is  treated  as  a  comment  line  in  the  input  file  and  is 
ignored  by  the  program.  Blank  lines  are  also  ignored  by  the  program. 


simulate  duration  npts  seed 

The  keyword  simulate  indicates  to  the  program  that  the  following  three  arguments 
in  sequence  are: 

•  duration:  Total  desired  duration  (in  seconds)  of  each  of  the  simulated 
wave  histories  (a  real  number). 


2.3.  Input  SpeciScation 
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•  npts:  Number  of  points  required  in  each  of  the  simulated  wave  histories 
(an  integer). 

•  seed:  A  real  number  (131071.0  for  example)  for  generation  of  random 
numbers.  The  seed  may  be  changed  by  the  user  in  order  to  generate  a 
different  set  of  wave  histories.  The  number  should  be  between  1.0  and 
2^1. 

The  resulting  dt  (time  step)  in  the  simulated  histories  is  duration/ npts.  The  mxgrd 
variable  in  the  driver  program  controls  the  maximum  number  of  points  allowed  (Max¬ 
imum  npts  =  2xmxgrd)  in  a  simulation.  The  released  driver  program  has  mxgrd  = 
4096,  so  that  specified  npts  can  to  up  to  8192.  The  user  may  increase  or  decrease 
mxgrd  to  suit  his/her  needs. 


depth  value 

This  line  specifies  the  water  depth  at  the  site  of  interest.  The  keyword  is  depth  and 
value  is  a  real  number  indicating  the  water  depth.  The  units  (meters,  feet,  etc.)  of 
this  value  should  be  consistent  with  the  units  of  other  input  parameters. 


psd  psdtype  psdLparameters 

This  line  specifies  the  spectrum  type  to  be  used.  The  keyword  is  psd  followed  by 
its  arguments,  psdtype  may  be  one  of  the  following:  jonswap,  bimodal,  or  boxcar. 
If  any  other  word  is  specified  for  psdtype,  then  it  indicates  to  the  program  that  an 
input  spectrum  is  specified  in  a  file  whose  name  is  same  as  the  word  specified  in  place 
of  psdtype.  More  details  regarding  this  input  specification  are  given  in  the  following 
subsection. 


define  varlimit  value  [optional  command] 

If  included,  this  line  defines  a  constant  varlimit  whose  value  is  a  real  number  (be¬ 
tween  0.0  and  1.0)  equal  to  value.  A  warning  is  issued  by  the  simulation  routine,  if 
the  estimated  second-order  power  above  Nyquist  frequency  is  more  than  value  times 
the  first-order  power.  If  this  line  is  not  provided  in  the  input  file,  then  a  default  value 
of  0.01  is  assigned  to  varlimit. 


define  gravity  value  [optional  command] 

If  included,  this  line  specifies  the  acceleration  due  to  gravity  in  consistent  units. 
value,  a  real  number,  is  assigned  to  gravity.  If  not  included,  a  default  value  of  9.807 
meters/sec^  is  assumed. 
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dciiiis  oingmsx  value  [Options,!  Commsnd]  If  includod,  this  lin6  spocifiss  tho  ni3.x- 
imum  frequency  in  radians  per  second  which  will  be  included  in  the  second  order 
transfer  function  calculation.  Specifically,  if  {cvm  +  ^n)  in  equation  2.5  is  greater  than 
omgmax,  then  that  (m,  n)  contribution  to  is  not  included  in  the  sum  in  equation 
numeq  x2. 

write  history  filenamel  filename2  [optional  command] 

If  included,  this  line  specifies  the  files  which  contain  the  simulated  time  histories  at 
each  location.  The  keyword  is  write  history.  The  file  filenamel  contains  the  un¬ 
derlying  first-order  (Gaussian)  histories,  and  filename2  contains  the  corresponding 
second-order  wave  histories.  The  format  of  the  output  is  presented  in  the  following 
section.  If  this  line  is  not  provided  in  the  input  file  then  default  names  ofgauss.hist 
and  ngauss.hist  are  assigned  to  the  output  first-  and  second-order  histories,  respec¬ 
tively. 


write  statistics  filenames  filename^  [optional  command] 

This  line  specifies  the  files  to  which  the  statistics  (mean,  standard  deviation,  skew¬ 
ness,  kurtosis,  minimum,  and  maximum  value)  of  the  simulated  histories  should  be 
written.  The  keyword  is  write  statistics.  The  statistics  for  the  simulated  first-order 
histories  at  each  spatial  location  is  written  out  in  filenames  and  the  statistics  for  the 
total  second-order  histories  are  written  in  filenamei-  If  this  line  is  not  provided  in 
the  input  file  then  default  names  of  gauss-stat  and  ngauss.stat  are  assigned  to  the 
output  files  for  the  first-  and  second-order  history  statistics,  respectively. 


location  nloc 

This  line  specifies  the  number  of  spatial  locations  at  which  both  the  first-  and  total 
second-order  wave  histories  should  be  simulated  in  time.  The  keyword  is  location. 
nloc  (an  integer)  specifies  the  number  of  locations  (maximum  location  allowed  is  50). 
valuel,  value2,  . . . ,  valuenloc  are  the  spatial  values  (real  numbers)  in  consistent  units. 
The  number  of  values  should  be  equal  to  nloc.  The  user  is  forewarned  that  the  spec¬ 
ification  of  the  spatial  locations  should  be  at  the  end  of  all  other  inputs  required. 


2.3.1  Wave  Spectrum  Specification 


The  input  wave  spectrum  can  be  specified  in  various  ways:  spectrum  values  in  an 
input  file,  or  spectrum  type  with  related  parameters  from  a  default  library. 


2.4.  Output  Format 
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psd  filename 

This  specifies  that  the  input  spectrum  be  read  from  an  input  file  named  filename.  The 
format  of  this  file  is  two  free-formatted  values  per  line,  specifying  a  natural  frequency 
in  radians  per  second  and  the  one-sided  PSD  (power  spectral  density)  ordinate  for 
that  frequency  in  consistent  units  of  squared  amplitude  (feet^,  meters^,  etc.)  per 
rad/sec.  Lines  beginning  with  a  “#”  are  regarded  as  comments  and  ignored.  The 
input  spectrum  may  be  specified  on  an  irregular  mesh.  This  is  internally  linearly 
interpolated  to  a  spectrum  on  a  regular  mesh  specified  by  duration  and  npts.  The 
spectral  ordinates  below  the  minimum  frequency  and  above  the  maximum  frequency 
specified  in  the  irregular  mesh  are  assumed  to  be  zero. 

Spectral  models  from  the  library  are  called  upon  using  any  one  of  the  following 
reserved  names  followed  by  their  parameters  (that  are  real  numbers): 


psd  jonswap  Hg  Tp  j 
psd  bimodal  Hg  Tp 
psd  boxcar  wio 


The  keyword  jonswap  invokes  a  JONSWAP  spectrum  parameterized  by  the  sig¬ 
nificant  wave  height  Hs  (defined  as  four  times  the  standard  deviation  of  the  wave 
elevation  process),  spectral  peak  period  Tp  (in  seconds),  and  the  peakedness  factor  7. 

The  keyword  bimodal  invokes  a  spectral  model  proposed  by  Torsethaugen  (Bitner- 
Gregerson  and  Haver,  1991).  This  subdivides  the  Hg-Tp  scattergram  into  three  re¬ 
gions,  and  assigns  bimodal  spectral  shapes  in  several  of  these  regions.  Therefore,  the 
only  input  required  for  this  bimodal  option  is  Hg  and  Tp. 

Finally,  the  keyword  boxcar  invokes  a  simple  band-limited  white-noise  model  of 
the  first-order  wave  spectrum.  Its  parameters  are  the  rms  a,,,  lower  cutoff  frequency 
a^io,  and  upper  cutoff  frequency  a;hj  of  the  first-order  wave  spectrum.  As  in  other  cases 
(e.g.,  the  user-defined  spectrum  at  various  frequencies),  the  frequencies  u\o,  and  ww 
here  are  assumed  to  be  in  units  of  rad/sec.  Note  also  that  non-zero  values  of  the  PSD 
at  zero  frequency  are  not  allowed  in  either  file  input  or  library  model  selections. 


2.4  Output  Format 


A  total  of  four  output  files  are  produced  by  the  driver  program.  Two  output  files 
contain  the  time  histories;  one  for  the  underlying  first-order  wave  histories  and  the 
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other  for  the  total  second-order  histories.  The  other  two  output  files  contain  wave 
statistics:  first  four  moments,  minimum  and  maximum.  Again,  results  for  the  first- 
and  second-order  wave  histories  are  separated  into  two  files. 


2.4.1  Time  History  Output 


As  noted  in  the  previous  section,  by  default  the  first-  and  second-order  histories  are 
written  to  the  files  gauss.hist  and  ngauss.hist.  Other  choices  of  output  filenames 
can  be  specified  by  the  optional  write  history  command.  The  format  of  this  output 
depends  on  the  number  of  spatial  locations  specified.  If  the  number  of  locations  (nloc) 
is  less  than  or  equal  to  8  then  the  output  is  in  Format  1  otherwise  the  output  is  in 
Format2.  Both  of  these  formats  write  out  3  header  lines  beginning  with  a  sign. 
These  are  to  be  treated  as  comment  lines  in  the  output  file.  In  order  to  explain  the 
two  format  styles,  say  that  the  spatial  locations  specified  are  xl,  x2,  x3,  . . . ,  xnloc. 

Formatl  outputs  data  in  n/oc-l-1  columns.  The  length  of  each  column  is  equal 
to  the  number  of  points  desired  in  each  simulation.  The  first  column  contains  the 
time  increments  in  seconds  going  from  0  to  T  with  dt  =  T/npts.  Columns  2  through 
nloc-\-l  contain  the  simulated  history  values  at  the  specified  locations  xl,  x2,  x3, 
xnloc,  respectively.  Thus,  column  2  contains  the  wave  elevation  at  location  xl, 
column  3  contains  wave  elevation  at  location  x2,  and  so  on. 

Format2  is  for  handling  nloc  greater  than  8.  The  output  begins  with  the  time 
increment  Ti  in  seconds  on  a  line  by  itself.  The  time  history  values  for  the  specified 
spatial  locations  at  time  Ti  are  written  in  the  next  line  onwards,  in  sets  of  10.  So  if  9 
locations  were  specified  (i.e.,  nloc  —  9)  then  the  time  increment  is  printed  on  a  line  by 
itself  followed  by  a  line  containing  9  time  history  values  at  that  time  increment.  The 
next  line  contains  the  next  time  increment  followed  by  another  set  of  9  values,  and 
so  on.  If,  on  the  other  hand  say  28  locations  were  specified,  then  a  time  increment  is 
written  on  a  line  followed  by  28  time  history  values  (corresponding  to  28  locations  at 
that  time  increment)  in  the  next  3  lines.  The  first  line  of  the  3  lines  contains  10  time 
history  values  for  the  first  10  locations  specified.  The  next  line  contains  10  history 
values  for  locations  11  through  20  and  the  following  line  which  is  the  third  line  of  the 
set  will  contain  only  8  history  values  for  location  21  through  28. 


2.5.  Example 
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2.4.2  Wave  Statistics  Output 


The  statistics  of  the  simulated  histories  are  also  estimated  by  the  driver  program. 
These  statistics  include  the  mean,  standard  deviation,  skewness,  kurtosis,  minimum, 
and  maximum.  As  noted  in  the  previous  section,  first-  and  second-order  simulation 
results  are  written  by  default  to  the  files  gauss.stat  and  ngauss.stat,  respectively. 
The  optional  command  write  statistics  can  alter  this  choice  of  output  filenames. 

The  output  format  in  both  of  these  files  begins  with  2  header  lines,  each  of  which 
begins  with  a  sign.  The  output  is  in  seven  columns.  The  first  column  specifies 
the  spatial  location.  The  following  six  columns  contain  statistics  of  the  wave  history 
at  the  spatial  location  specified  in  column  1.  Columns  2  through  7  contain,  the  mean, 
standard  deviation,  skewness,  kurtosis,  minimum, ,  and  maximum  value  in  that  order. 

The  next  section  presents  some  sample  output  files,  to  illustrate  use  of  the  WAVEMAKER 
routine. 


2.5  Example 


In  this  section,  we  present  a  sample  problem  (copies  of  input  and  output  files  are 
enclosed  on  disk).  To  illustrate,  consider  a  simulation  which  samples  the  wave  process 
at  regular  intervals  of  length  dt=0.5  [sec]  over  a  total  duration  of  T=2048  [sec];  i.e., 
npts=4096  points.  We  assume  here  the  first-order  wave  spectrum  to  be  of  JONSWAP 
form,  with  Hg  =  12  [m],  7^  =  14  [sec]  and  a  peakedness  faotor  7  =  3.3.  We  further 
seek  to  generate  wave  histories  at  2  spatial  locations:  0  and  60  [m].  Our  input  file  for 
simulating  waves  using  WAVEMAKER  should  be: 
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#  Gaussian  and  Nongaussian  Wave  Input  File 

simulate  2048.0  4096  8123872.0 
depth  70.0 

psd  jonswap  12.0  14.0  3.3 

write  history  gauss.hist  ngauss.hist 

write  statistics  gauss.stat  ngauss.stat 

define  varlimit  0.01 

define  gravity  9.807 

define  omgmax  4.2 

location  2 

0.0 

60.0  _ 


Alternatively,  if  we  intend  to  use  the  default  definitions  in  the  program  then  our 
input  file  could  be  (this  will  produce  the  same  output  as  the  extended  version  of  the 
input  file): 


#  Gaussian  and  Nongaussian  Wave  Input  File 

#  (Alternative  format) 

simulate  2048.0  4096  8123872.0 
depth  70.0 

psd  jonswap  12.0  14.0  3.3 
define  omgmax  4.2 
location  2 
0.0 

60.0  _ 


The  output  files  created  are:  gauss.hist,  ngauss.hist,  gauss.stat,  and  ngauss.stat. 
The  contents  of  these  are  listed  in  the  following  table: 


Output  File 

Contents 

gauss.hist 

ngauss.hist 

First-Order  Time  History 
Second-Order  Time  History 

gauss.stat 

ngauss.stat 

First-Order  History  Statistics 
Second-Order  History  Statistics 

Portions  of  each  of  these  output  files  are  given  in  the  Appendices. 


2.5.  Example 
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Figures  2.2  and  2.3  show  a  comparison  of  the  simulated  histories.  Figure  2.2  com¬ 
pares  the  simulated  first-  and  second-order  wave  time  histories  at  the  same  location 
(the  first  of  the  two  requested,  defined  arbitrarily  as  a:=0  [m]).  The  file  ngauss.stat 
includes  the  estimated  skewness  of  the  second-order  waves.  At  this  location  the 
second-order  wave  history  has  a  skewness  of  about  0.2.  This  positive  skewness  (com¬ 
pared  to  zero  skewness  of  Gaussian  waves)  indicates  the  systematic  nonlinear  eflfects. 
This  also  gives  rise  to  an  asymmetry  between  peaks  and  troughs;  in  particular,  the 
extreme  wave  crest  in  the  second-order  simulation  systematically  exceeds  the  corre¬ 
sponding  extreme  trough  in  absolute  value.  This  tendency  may  be  significant  for 
potential  deck  impact  problems,  particularly  in  older  jacket  structures  with  relatively 
low  deck  levels. 

Figure  2.3  compares  the  second-order  wave  histories  predicted  by  simulation  at 
the  two  spatial  locations,  separated  by  60  [m].  The  phase  lag  at  these  two  locations 
is  evident  in  the  plot.  It  can  be  seen  that  a  crest  height  at  0  [m]  does  not  necessarily 
imply  a  crest  of  the  same  height  at  60  [m]. 


Elevation  (meters)  ^  Elevation  (motors) 
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'igure  2.2:  Simulated  first-  and  second-order  wave  histories  at  location  0.0 


Time  (sec) 

Figure  2.3:  Simulated  second-order  wave  histories  at  location  0.0  and  60.0 


Chapter  3 


Identification  of  First-Order  Waves 


3.1  Introduction 


In  ocean  engineering  practice  it  is  common  to  assume  the  waves  to  be  Gaussian 
and  any  nonlinearity  in  the  waves  is  embedded  in  the  structural  response  analysis 
(e.g.,  WAMIT,  1995).  It  has  been  shown  in  Winterstein  and  Jha,  1995  that  observed 
time  histories  generally  contain  nonlinearities,  it  is  thus  imperative  to  remove  any 
second-order  effects  in  the  incident  waves  so  that  we  do  not  double-count  these  effects 
in  the  resulting  response  estimation.  Recent  studies  (Ude  and  Winterstein,  1996) 
have  demonstrated  the  impact  of  double-counting  such  second-order  effects  on  various 
structural  response  characteristics. 

The  methodology  to  identify  the  underlying  first-order  waves  is  to  seek  the  implied 
first-order  wave  history  which,  when  run  through  the  second-order  wave  predictor, 
yields  an  incident  wave  that  agrees  with  the  target  observed  history  at  each  time 
point.  This  identification  is  performed  using  a  Newton-Raphson  scheme  to  achieve 
simultaneous  convergence  at  each  complex  Fourier  component.  If  the  observed  history 
has  N  components,  we  iteratively  solve  N  simultaneous  nonlinear  equations  to  identify 
the  first-order  components. 

Due  to  computer  memory  limitations,  the  identification  of  the  first-order  history 
is  performed  on  short  contiguous  windows  of  the  observed  history.  This  window  size 
(mien)  can  be  made  equal  to  the  observed  history  length  in  WAVEMAKER  if  the  computer 
has  sufficient  RAM  and  swap  space. 
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3.2  Methodology 

The  idea  here  is  to  identify  the  implied  first-order  history  T]i{t)  (of  an  observed  history 
Vobsit))  which,  when  run  through  the  second-order  predictor,  yields  an  incident  wave 
that  agrees  with  Tjobs{t).  In  the  first-order  wave  process  771  (t)  (see  Eq.  2.7),  written  as 
a  Fourier  sum  of  N  frequencies, 

N/2  N 

Vli^)  COs{uikt  -I-  9k)  =  ^2  (3-1) 

k=l  k=l 

we  need  to  identify  only  the  lower  half  Xk  components,  since  the  upper  half  values 
are  complex  conjugates  of  the  lower  half.  Let  us  denote  Xk  =  Uk  +  iVk,  where  Uk,  14 
are  the  real  amd  imaginary  parts  of  the  complex  Fourier  component  Xk,  respectively. 

The  predicted  second-order  wave  process  (see  Eq.  2.11)  as  evaluated  from  the 
QTFs  is 

N/l  N/2 

Am(t)  =  2Re  E  E  (3.2) 

T7i=l  n=l 

This  may  be  rewritten  in  the  form  of  a  Fourier  sum  as 

A772(t)  =  El4e'‘"'‘‘  (3-3) 

ifc=i 

where  Yk  =  -I-  Fyf  are  the  combined  sum  and  difference  frequency  components. 

Here  again,  I4  possesses  conjugate  symmetry  so  that  only  the  lower  half  contains 
unique  information.  Y/^  can  be  shown  to  be 

=  E 

m+n,/: 

=  E  -  K,V„)  +  (3.4) 

m+n,A: 

where  the  summation  symbol  indicates  a  double  summation 

N/2  N/2 

^  ^  such  that  Um  +<x)n=^k  (3-5) 

m-hfiik  m=l  n=l 

and 

n-  =  E 

m’-n,k 

=  E  mmUn+VM+iiVmUn-UmVn)]H-, 

m—riyk 


(3.6) 
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where 


N/2  N/2 

^  £  such  that  \u}m  - 

m—n,k  771=1  n=l 


The  combined  predicted  wave  process  is 


Vpredit)  —  Vli.l')  *i" 


(3.7) 


(3.8) 


The  identification  scheme  strives  to  simultaneously  match  VedC*)  observed 

wave  history  7;obs(i)  at  every  value  of  t.  Alternatively,  we  can  perform  the  identi¬ 
fication  in  the  frequency  domain  and  strive  to  simultaneously  match  the  predicted 
Fourier  components  to  the  observed  Fourier  components  at  all  frequencies. 


Voheif)  can  be  represented  in  the  frequency  domain  as 

T;obs(t)  =  E  (3.9) 

*=1 

where  ZkS  also  possess  conjugate  symmetry.  If  the  first-order  components  are  iden¬ 
tified  exactly,  from  Eq.s  3.1,  3.3  and  3.9  we  will  have 

Zk  =  Xk  +  Yk  ;  foralljk  =  l...iV/2  (3.10) 

Note  that  the  upper  half  values  can  be  obtained  from  conjugate  symmetry  of  the 
lower  half  values.  In  the  Newton-Raphson  identification  scheme  we  will  try  to  simul¬ 
taneously  minimize  Xk  +  Yk  -  Zk;  for  k  =  1 ...  N/2  to  achieve  convergence.  Now,  this 
scheme  requires  a  Jacobian  of  Xk  +  Yk-  Zk  with  respect  to  the  unknowns  X^-such  a 
complex  differentiation  will  lead  to  numerical  discontinuities  so  we  will  minimize  an 

equivalent  real  function  \/Ef  fk/^  instead,  where  for  fc  =  1 . . .  N/2 

fk  =  Re{Xk  +  Yk-Zk) 

fk+N/2  =  lm{Xk  +  Yk-Zk)  (3.11) 


The  identification  of  the  lower  half  Xk  values  requires  a  simultaneous  solution  of 

the  nonlinear  equations  in  3.11  such  that  fk->0  for  all  A:  =  1 . N,  or  alternately 

fk  /X  — >  0.  We  will  formulate  the  Newton-Raphson  scheme  in  vector  form  as 


■  ReX  ■ 

-H 

■  ReY  ■ 

ReZ 

ImX 

ImY 

ImZ 

(3.12) 


where  bold  face  letters  denote  vectors,  and  vectors  X,  Y,  Z  contain  the  complex 
Fourier  components  Xk,  Yk,Zk,k-\...  N/2,  respectively.  Here,  [-^-]  is  a  vector 
containing  the  real  part  of  X  in  the  upper  half  and  the  imaginary  part  of  X  in  the 
lower  half. 
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Let  us  denote 


u  ■ 

V 


B  = 

C  = 

Note  that  the  vector  A,  of  length  N,  is  constructed  such  that  lower  half  values 
are  the  real  parts  of  X*;  A:  =  1 . . .  N/2  and  the  upper  half  is  the  imaginary  part 
of  X*;  k  =  I...N/2.  Similarly,  B  and  C,  each  of  length  N,  contain  real  and 
imaginary  parts  of  the  lower  half  of  the  second-order  correction  and  the  observed 
Fourier  components,  respectively.  The  elements  of  A  and  B  are  denoted  by  a/  and 
bk,  respectively,  where  l,k  =  The  objective  function  in  vector  notation  now 

f(A)  =  A-\-B-C  (3.14) 

A  first-order  Taylor  approximation  of  f(A)  about  a  given  A(°)  is 


ReY 


ImY 

ReZ 

ImZ 


(3.13) 


■  ReX  ■ 

ImX 

f  (A)  =  f  (a^®))  +  [J]  (a  -  AW)  (3.15) 


where  [J]  is  a  NxN  Jacobian  matrix  denoting  the  derivatives  of  the  elements  fk  in 
vector  f(A)  with  respect  to  each  of  the  unknowns  a/  in  A  where  k,l  =  The 

Newton-Raphson  scheme  at  iteration  p  -M  is  then  formulated  as 

A(p+i)  =  a(p>  -F  h  (3.16) 


where  h,  a  vector  of  length  N,  is  found  from  a  Cholesky  decomposition  followed  by 
a  back-substitution  scheme  from 

[J]h  =  -f  (a(p))  (3.17) 

It  can  be  easily  shown  from  Eq.  3.14  that  the  entries  Jk,i  of  the  matrix  [J]  are 

where  dbk/dai  indicates  the  partial  derivative  of  6*  with  respect  to  ai,  and 

5ki  =  [^  (3.19) 

(0  otherwise 

To  find  dbk/dai,  recall  from  notation  in  3.13 

bk  =  Rel*  and  bk+N/2  —  fo'^  A:  =  1 . . .  N/2 

ai  =  ImX;  =  Ui  and  ai+N/2  =  ImX;  =  Vi  for  /  =  1 . . .  N/2 
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so  that  from  Eq.s  3.4  and  3.6  we  have 

aReFfe 


dUi 


m+n,fc  m—Tijk 


=  E  +  E  (Ki™i  +  v„WH»n  (3.20) 

aVi  ^ 

51ml*; 
dUi 


TO+n,/t  m-n,A; 

=  53  +  K»^mi)  -^mn  "I"  (Im^n/  “  ^n^ml)  H„ 


m-^Tiyk 


m—riyk 


=  £  (P„5™,  +  £4.3.^)  HI.  +  E 

oVi 


m+Jiyk 


T7l““Tl}A/ 


Schematically, 


[J]  =  [/]  + 

where  [/]  is  the  identity  matrix. 


- 1 

dReYk  1 
dV, 

5Imljfe 

dUi 

Q: 

1 _ 

(3.21) 


3.2.1  Ramp 


In  general,  the  observed  input  wave  time  series  can  start  or  end  from  any  value.  If 
the  two  ends  of  the  wave  record  or  of  the  window  being  processed  are  not  continuous 
in  both  offset  and  slope,  the  end  discontinuity  will  introduce  high  frequency  ’’noise” 
in  the  predicted  time  series,  as  is  further  discussed  in  chapter  4.  The  discontinuities, 
and  so  the  high  frequency  ”  noise,”  can  be  suppressed  by  prepending  and  appending 
a  ramping  sequence  to  the  wave  history  in  each  window. 

The  ramp  is  linearly  raised  from  or  collapsed  to  the  mean  of  the  input  using  a 
modulated  sinusoidal  function: 

R{t  +  U)=x  +  cos{u}t)  +  —  sin(a;t)]  {-td<t<0)  (3.22) 

td  ^ 

in  which  x  is  the  mean  of  the  observed  input  time  series  and  Xq  =  a;(0)  —  x  and 
Vq  =  (x(l)  -  x{0))/dt  -  xo/td  are  the  initial  value  and  modified  velocity  of  the  zero 
mean  input  x.  dt  is  the  sampling  time  interval.  The  frequency  oj  of  the  ramp  is 
estimated  as  the  average  frequency  of  the  first  and  last  ten  cycles  of  the  input  time 
series.  Once  the  required  output  time  series  has  been  identified  with  the  augmented 
input  time  series,  its  initial  and  final  parts  corresponding  to  the  ramps  are  removed. 
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3.2.2  Newton-Raphson  Scheme 

The  algorithm  for  the  Newton-Raphson  scheme  followed  in  WAVEMAKER  is 


1.  Estimate  C  from  observed  history  (Eq.s  3.9,  3.13) 

2.  Initial  Guess  A  =  C 

3.  Estimate  B  from  A  (Eq.s  3.4,  3.6,  3.13) 

4.  Find  f(A)  (Eq.  3.14) 

5.  Find  [J]  (Eq.s  3.20,  3.21) 

6.  Solve  [.7]h  =  -f(A)  to  find  h 

7.  Update  A  =  A  -f  h 

8.  Check  Convergence  (see  next  section): 

If  converged  terminate  else  go  to  3 


3.2.3  Convergence  Criteria 


The  Newton-Raphson  iteration  scheme  is  terminated  based  on  the  following  condi¬ 
tions: 


•  Program  Converged:  If  the  rms  of  the  increment  vector  h  =  ^l/N 
is  less  than  a  specified  tolerance,  the  program  is  said  to  have  converged. 
This  convergence  tolerance  is  specified  as  a  fraction  a  (=  0.0001  in  WAVEMAKER) 
of  the  standard  deviation  of  the  observed  wave  history  cr,,,obs- 

•  Program  Diverging:  If  the  rms  of  the  identified  first-order  history 

at  any  iteration  p  is  larger  than  a  specified  fraction  jS  (=  200  in  WAVEMAKER) 
of  o-„,obs  then  the  identification  scheme  is  restarted  with  a  smaller  initial 
guess  which  is  a  truncated  and  scaled  down  version  of  C.  The  truncation 
point  is  at  twice  the  peak  spectral  frequency  of  77obs(t)  and  the  scaling 
factor  is  factnu^  (=  0.9’’  in  WAVEMAKER),  where  r  is  the  number  of  restarts 
needed  so  far.  Thus  the  restart  guess  in  complex  Fourier  notation  is 


Xk 


facttltl  Zh  ^  UJk  ^  ^a^peak 
0  ;  otherwise 


(3.23) 


•  Maximum  Iterations  Reached:  If  the  maximum  allowed  iterations, 
specified  by  the  variable  mxiter  (=  10  in  WAVEMAKER),  is  reached  and  the 
program  has  still  not  converged,  then  the  program  restarts  the  Newton- 
Raphson  identification  scheme  with  a  smaller  initial  guess  =  factnvfC. 

•  Maximum  Restarts  Reached:  If  the  maximum  allowed  restarts,  speci¬ 
fied  by  the  variable  nuiter  (=  5  in  WAVEMAKER),  is  reached  then  the  program 
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terminates  the  identification  scheme  in  the  present  window  and  proceeds 
to  identify  in  the  next  observed  history  window. 


3.2.4  Implementation 

The  first-order  components  for  the  observed  wave  history,  of  length  A^obs>  are  identified 
in  contiguous  windows,  each  of  length  N  <  Nobs-  The  identification  analysis  is 
performed  in  this  way  to  minimize  the  computer  memory  usage  by  WAVEMAKER.  Recall 
that  the  Jacobian  matrix  [J]  is  a  NxN  matrix  and  the  memory  usage  is  directly 
governed  by  the  matrix  size  of  this  variable.  In  principle,  if  there  is  sufficient  memory 
we  could  set  N  =  Nobs  and  identify  the  first-order  component  for  the  entire  observed 
history  in  one  window,  however,  this  is  not  usually  not  the  case  and  we  resort  to 
identifying  in  contiguous  windows,  as  shown  in  Fig.  3.1. 

The  first-order  component  is  identified  independently  in  each  of  the  windows  in 
sequence.  The  last  window  is  skipped  if  its  contains  points  less  than  N. 


3.3  Input  Specification 


The  input  specification  for  the  identification  of  first-order  wave  process  is  in  a  command¬ 
line  format  similar  to  the  simulation  input.  A  typical  input  file  for  identification  is: 


#  Typical  input  file:  syntax  description 

identify  filename  dt  winsize  nramp 

depth  value 

define  varlimit  value 

define  gravity  value 

define  omgmax  value 

write  history  filenamel  filename2 


#  Typical  input  file:  syntax  description 

Any  line  beginning  with  a  symbol  is  treated  as  a  comment  line  and  is  ignored. 
Blank  lines  in  the  input  file  are  ignored,  as  well. 

define  varlimit  value 
define  gravity  value 
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Figure  3.1:  Identification  of  first-order  wave  components  is  done  in  contiguous  win¬ 
dows  of  the  observed  history 


The  keywords  varlimit,  gravity  and  omgmax  have  the  same  meaning  as  in  the 
simulation  section  and  the  user  is  referred  to  this  section  to  understand  the  usage  of 
these  commands. 

depth  value 

The  keyword  depth  as  in  the  simulation  section  indicates  the  water  depth  at  which 
the  identification  analysis  is  to  be  performed. 

identify  filename  dt  winsize  nramp 

The  keyword  identify  indicates  to  the  program  that  the  user  intends  to  identify  the 
underlying  first-order  wave  history  for  a  given  observed  wave  history.  This  command 
requires  three  arguments  which  in  sequence  are: 

•  filename:  The  name  of  the  file,  a  character  string,  containing  the  observed 
wave  time  history  for  which  the  underlying  first-order  wave  history  is  to 
be  identified.  The  data  in  the  first  column  in  filename  is  read  as  the 
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observed  wave  time  history.  Any  blank  lines  in  filename  or  lines  that  do 
not  begin  with  a  number  are  ignored. 

•  dt:  This  value,  a  real  number,  indicates  the  time  resolution  of  the  wave 
history  provided  in  filename.  In  other  words,  dt  is  the  time  difference 
between  two  successive  elevation  values  in  the  observed  wave  history. 

•  winsize:  An  even  integer  value  indicating  the  window  size  or  the  number 
of  points  of  the  provided  wave  history  to  be  used  in  each  Newton-Raphson 
iteration.  The  first-order  wave  components  are  identified  in  windows  (of 
size  winsize)  in  sequence  for  the  provided  time  history.  If  the  last  win¬ 
dow  contains  number  of  points  less  than  winsize,  then  winsize  is  reset  to 
the  remaining  number  of  points.  If  the  original  history  contains  an  odd 
number  of  points  then  the  last  point  is  ignored. 

•  nramp:  An  integer  value  indicating  the  number  of  points  to  be  included  at 
each  end  of  the  portion  of  the  wave  history  in  each  window.  The  number 
of  data  points  processed  in  each  window  will  be  mnsize—2xnramp.  If 
nramp  is  not  included  after  the  keyword  identify  or  is  set  to  zero,  then  no 
ramps  are  included. 

The  maximum  value  of  winsize  is  mien  set  to  512  points  in  WAVEMAKER  and  can 
be  changed  according  to  the  user’s  needs  or  according  to  the  computer’s  limitations. 
Note  that  we  require  mien  <  2xmxgrd  in  the  program.  These  dimension  values  are 
set  in  this  way  so  as  to  minimize  the  memory  requirements  of  WAVEMAKER. 

•write  history  filenamel  filenameB 

The  command  write  history  is  used  to  specify  the  file  names  where  the  identi¬ 
fied  histories  are  to  be  written.  The  identified  first-order  wave  history  is  written  in 
file  filenamel  and  the  identified  second-order,  combined  first-  and  second-order,  and 
the  observed  wave  histories  are  written  in  file  filename2.  Default  values  assigned  to 
filenamel  and  filenameB  are  gauss.hist  and  ngauss.hist,  respectively. 


3.4  Output  Format 


The  output  file  names  are  governed  by  the  command 
•write  history  filenamel  filename2 

with  default  names  being  gauss.hist  and  ngauss.hist  for  filenamel  and  filename2, 
respectively. 


26 


Chapter  3.  Identification  of  First-Order  Waves 


filenamel  contains  the  identified  first-order  wave  history  in  a  column  of  real  num¬ 
bers.  Each  line  of  this  file  contains  one  real  number  indicating  the  elevation  of  the 
first-order  wave  history  (see  example  output  files  in  Appendix  B).  The  time  resolu¬ 
tion  of  this  first-order  history  is  dt,  equal  to  the  dt  provided  in  the  input  file  using  the 
command  identify.  This  file  also  contains  comment  lines  that  begin  with  a  “#”  sym¬ 
bol  as  the  first  character  on  the  line.  The  first  comment  line  contains  information  on 
the  contents  of  the  file,  and  the  following  comment  lines  contain  3  integers:  the  first  is 
the  window  number  being  identified,  the  second  is  the  number  of  iterations  required 
for  convergence,  and  the  third  is  the  number  of  restarts  needed  for  convergence. 

filename2  contains  the  second-order  correction,  the  combined  first-  and  second- 
order  waves,  and  the  observed  wave  time  history.  The  second-order  correction  is  found 
from  the  identified  first-order  waves,  and  these  two  are  added  together  to  yield  the 
combined  second-order  wave  history.  These  histories  are  provided  in  three  columns 
in  filename2,  or  in  other  words  each  line  contains  three  real  numbers:  the  first  is  the 
second-order  wave  elevation,  the  second  is  the  combined  identified  wave  elevation, 
and  the  third  is  the  observed  wave  history  (see  example  output  file  in  Appendix  B). 
A  match  of  the  total  identified  and  the  observed  wave  histories  will  verify  successfuly 
identification  by  WAVEMAKER.  The  time  resolution  of  each  of  these  histories  is  dt. 
This  file  also  contains  comment  lines  beginning  with  a  “#”  symbol  that  provides 
information  similar  to  the  comment  lines  in  filenamel. 


3.5  Examples 


In  this  section  we  present  two  sample  problems  to  illustrate  the  use  of  the  identifi¬ 
cation  capabilities  of  WAVEMAKER.  Example  1  is  based  on  the  example  presented  in 
the  simulation  chapter.  Sample  input  and  output  files  of  this  identification  example 
are  included  in  the  distribution  diskette.  Example  2  presented  here  demonstrates 
the  identification  of  first-order  components  of  a  measured  wave  tank  history.  Note 
that  sample  input  or  output  files  of  this  second  example  are  not  included  in  the 
distribution. 


3.5.1  Example  1 


The  example  of  the  simulation  capabilities  of  WAVEMAKER  involved  simulating  a  second- 
order  wave  history  characterized  by  a  JONSWAP  spectrum  with  Hs  =  12  [m],  Tp  = 
14s  and  7  =  3.3  in  70  [m]  water  depth.  We  will  use  the  combined  second-order 
simulated  history  and  try  to  identify  its  first-order  wave  component  and  compare  it 
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to  the  input  first-order  component  used  to  simulate  the  combined  wave  history.  The 
input  file  for  the  identification  run  is 


#  Wave  Identification  Input  File 

identify  hist.dat  0.5  512 
depth  70.0 

write  history  gauss.ide  ngauss.ide 
define  varlimit  0.01 
define  gravity  9.807 
define  omgmax  4.2 _ 


Alternatively,  if  we  intend  to  use  the  default  definitions  in  the  program  then  our 
input  file  could  be  (this  will  produce  the  same  output  as  the  extended  version  of  the 
input  file): 


#  Wave  Identification  Input  File 

#  (Alternative  format) 

identify  hist.dat  0.5  512 

depth  70.0 

define  omgmax  4.2 

write  history  gauss.ide  ngauss.ide 


The  input  file  hist.dat  contains  a  column  of  real  numbers  (see  sample  files  listed 
in  the  appendix)  which  will  be  read  in  as  the  observed  wave  history.  The  first- 
order  components  will  be  identified  for  this  observed  history  and  placed  in  the  file 
gauss.ide.  The  corresponding  second-order  components,  the  combined  identified 
history  and  the  observed  wave  history  are  written  in  the  file  ngauss.ide. 


Output  File 

Contents 

gauss.ide 

ngauss.ide 

Identified  First-Order  Time  History 
Corresponding  Second-Order  Time  History 

Figure  3.2  shows  the  observed  wave  spectrum  and  the  identified  first-order  spec¬ 
trum  along  with  the  corresponding  second-order  wave  spectrum.  We  see  that  small 
second-order  contribution  to  the  power  spectrum,  roughly  a  decade  below  the  first- 
order  spectrum  even  at  frequencies  twice  the  peak  spectral  frequency,  suggests  the 
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difficulty  in  identifying  these  components.  Figure  3.3  shows  the  observed  wave  his¬ 
tory  and  the  identified  first-order  wave  history  in  cycles  around  the  maximum  crest 
height.  Compare  this  to  the  simulation  example  where  we  solve  the  forward  problem 
of  finding  the  combined  (first-  plus  second-order)  history  from  a  given  underlying 
first-order  wave  spectrum.  The  identified  first-order  component  in  Fig.  3.3  is  almost 
the  same  as  the  underlying  first-order  component  (denoted  Gaussian)  in  Fig.  2.2  and 
these  two  are  shown  together  in  Fig.  3.4.  Note  how  close  the  two  first-order  compo¬ 
nents  are,  and  any  numerical  differences  can  probably  be  further  reduced  by  using  a 
larger  window  size  (greater  than  512,  for  example)  in  the  identification  scheme. 


3.5.2  Example  2 


In  this  example  we  will  identify  the  underlying  first-order  wave  component  for  a 
measured  wave  tank  history  that  reflects  a  water  depth  of  175m.  For  this  example 
the  measured  history  is  located  in  file  wave.dat  and  has  a  dt  =  0.3354  seconds.  We 
will  use  windows  of  winsize  =  512  to  identify  the  first-order  components.  The  input 
to  WAVEMAKER  is: 


#  Wave  Identification  Input  File 

#  using  default  definitions 

identify  wave.dat  0.3354  512  100 
depth  175.0 

define  omgmax  4.2 _ 


Figure  3.5  shows  a  portion  where  the  maximum  crest  height  occurs  in  the  mea¬ 
sured  wave  tank  history.  The  figure  also  shows  the  identified  first-order  and  the 
corresponding  second-order  wave  histories.  Note  how  the  second-order  wave  compo¬ 
nent  affects  the  first-order  peaks,  amplifying  the  crests  and  moderating  the  troughs. 
Figure  3.6  shows  the  wave  spectra  for  the  measured  history  along  with  the  first-order 
and  the  second-order  spectra.  Again,  observe  that  the  second-order  energy  is  signifi¬ 
cantly  small  compared  to  the  first-order,  however,  phase  locking  of  the  first-  and  the 
second-  component  (Fig.  3.5)  leads  to  larger  crests  and  flatter  troughs. 


spectrum  (m^2/Hz] 


Ffiingjig 
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Figure  3.2:  Wave  spectrum:  observed  vs.  identified  first-  and  second-orde 


Figure  3.3:  Wave  history:  observed  vs.  identified  first-  and  second-order 


Elevation  [m] 
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Chapter  4 


Prediction  of  Second- Order  Random 
Waves 


4.1  Introduction 


Chapters  2  and  3  of  this  document  outline  the  Simulate  and  Identify  options  of 
wavemaker.  The  Predict  option  builds  on  the  same  theory  and  implementation  to 
provide  another  capability:  to  propagate  an  observed  input  time  series  through  a 
second-order  model  to  produce  the  first-  and  second-order  outputs. 

Chapter  2  describes  a  procedure  for  simulating  a  second-order  random  wave  pro¬ 
cess  corresponding  to  a  user-specified  wave  power  spectral  density.  These  wave  histo¬ 
ries  include  a  second  order  Stokes  correction  at  the  sum  and  difference  of  each  wave 
frequency  pair,  as  calculated  in  the  frequency  domain.  An  inverse  FFT  is  applied  to 
this  frequency  domain  representation  of  the  wave  process  to  produce  a  time-history. 
Time  histories  for  multiple  spatial  locations  are  generated  by  applying  the  linear  dis¬ 
persion  relation  in  the  frequency  domain  of  the  underlying  first  order  wave  process 
prior  to  inclusion  of  the  second-order  correction  at  each  spatial  location. 

Chapter  3  describes  a  procedure  for  identifying  an  underlying  first-order  process 
from  a  user  specified  wave  history.  This  identification,  an  inverse  feature  to  sim¬ 
ulation,  is  based  on  a  Newton-Raphson  scheme  to  solve  N  simultaneous  nonlinear 
equations  to  identify  the  first-order  waves  which,  when  run  through  the  second-order 
wave  predictor,  matches  the  observed  waves. 

This  chapter  describes  a  third  execution  option  for  the  Wavemaker  program.  The 
predict  option  uses  the  underlying  first-order  wave  history  generated  by  the  identify 
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option  at  one  location  as  input  to  generate  consistent  first  and  second  order  wave 
time  histories  at  other  locations  using  essentially  the  same  theory  as  the  simulate 
option.  Thus,  the  Predict  option  combines  the  utility  of  the  procedures  developed  in 
Chapters  2  and  3. 


4.2  Overview 


Measured  wave  input  data,  such  as  that  collected  with  a  wave  probe,  is  first  identified 
as  described  in  Chapter  3.  The  Identified  underlying  first-order  wave  process  is  used 
as  the  input.  An  FFT  procedure  generates  a  discrete  Fourier  sum  representation  of 
the  wave  history  in  the  frequency  domain.  This  frequency  domain  representation  of 
the  wave  process  is  equivalent  to  that  generated  in  Chapter  2,  except  that  the  sum 
now  corresponds  to  a  user-specified  wave  time  history,  rather  than  a  user-specified 
power  spectral  density.  Time  histories  for  multiple  spatial  locations  are  generated  by 
use  of  the  dispersion  relation  with  a  second  order  correction  as  described  in  Chapter 
2. 


The  following  flowchart  shows  the  flow  of  the  calculation  for  the  predict  option. 
The  user  specifies  an  input  wave  history,  which  will  generally  be  a  measured  wave. 
Identify,  as  described  in  Chapter  3,  is  used  to  separate  the  total  process  into  consistent 
first  and  second  order  processes.  The  linear  dispersion  relation  is  applied  to  the 
consistent  first  order  part  of  the  wave  history,  which  produces  an  equivalent  first 
order  history  at  user  specified  locations.  The  predict  option  then  applies  the  inverse 
of  the  Identify  process  such  that  the  consistent  second  order  part  of  the  wave  history 
is  recreated  at  the  new  locations.  It  is  the  use  of  the  same  “black  box  in  the 
separation  (identify)  and  recreation  (predict)  calculations  that  ensures  the  newly 
created  second  order  part  of  the  wave  histories  at  the  new  locations  are  consistent 
with  those  identified  at  the  original  input  wave  location. 

The  implementation  of  the  prediction  capability  differs  from  simulation  only  in  the 
first  step.  In  simulation,  complex  Fourier  amplitudes  of  the  input  were  generated  from 
the  PSD  model.  In  prediction,  they  are  generated  by  simply  FFT’ing  the  observed 
input  time  series.  Once  the  input  is  described  by  complex  Fourier  amplitudes,  the 
procedure  is  identical  to  that  of  simulation. 

This  tool  provides  an  easy  way  to  compare  model  test  results  with  real-world 
observations  or  to  compare  specific  model  test  results  such  as  airgap  measurements 
at  a  known  location  with  a  model  test  reference  wave  history. 


4.3.  Methodology 
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Figure  4.1:  Overview  of  Prediction  Calculation  Flow 


4.3  Methodology 

4.3.1  Underlying  Theory  and  Assumptions 

We  first  consider  the  first-order  wave  elevation,  at  a  specific  reference  location 
(say  a:=0).  For  either  frequency-domain  analysis  or  time-domain  simulation,  it  is 
convenient  to  write  T)i{t)  as  a  discrete  Fourier  sum  over  positive  frequencies  Uk: 

Vi{t)  =  cos{u}kt  +  ^k)  =  Re  ^  (4.1) 

*=1  jfc=i 


For  purposes  of  prediction,  the  lowest  frequency  interval  du  is  governed  by  the 
total  period  T  of  the  simulation,  while  the  highest  frequency  is  governed  by  the  time 
interval  between  data  points  in  the  input  wave  history: 
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The  second-order  wave  at  this  elevation,  T]2{t),  is  calculated  from  T]i(t)  as 

n2{t)  =  T]i{i)  +  ^V2{t)  (4.3) 

in  which  includes  second-order  corrections  at  sums  and  differences  of  all 

wave  frequencies: 


AihW  =  R*  E  E  +  JZ+^e'#""-*"-*'***-****!]  (4.4) 

m=l n=l 

4.3.2  Implementation 

The  prediction  method  requires  an  input  time  history  representing  the  water  surface 
elevation  measurements.  The  measurements  are  required  to  be  at  everily  spaced  time 
intervals  to  take  advantage  of  discrete  FFT  techniques. 

An  FFT  is  applied  to  the  time  history  to  calculate  consistent  values  of  Ak  and  9k 
in  Eq.  2.1. 

Eq.  2.1  is  then  rewritten  as 

nptsl2  npts 

Vi{t)  =  E  cos(a;*t  9k)  =  Re  £  X.e^**  (4.5) 

fc=i  *=i 


Here  the  Xk  are  complex  Fourier  coefficients.  The  lower  half  of  these  directly 
reflect  values  of  amplitude  Ak  and  phase  9k  consistent  with  the  user  specified  wave 
history  at  frequency  0Jk=^  ‘ 

Xk  =  ;  k  =  l...npts/2  (4.6) 

The  upper  half  are  the  complex  conjugates  of  the  lower  half  due  to  the  double  sided 
nature  of  the  FFT  in  equation  (etal)  (the  symbol  “*”)  of  the  lower  half: 

Xnpu-k  =  X*k;  k  =  l...npts/2 


(4.7) 
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This  reflects  that  unique  information  is  contained  only  the  lower-half  frequencies; 
indeed,  any  information  in  the  upper  half  frequencies  (above  the  Nyquist)  is  obscured 
by  aliasing. 


4.4  End-Point  Continuity  of  a  Wave  Record 

4.4.1  Predict  Option 


Actual  measured  waves  are  generally  not  continuous  in  slope  or  oflfeet  between  the 
ends  of  the  measured  wave  record.  Subsequently,  the  identified  first  order  wave  process 
used  as  the  input  to  the  Predict  option  is  also  discontinuous  between  its  ends.  This 
discontinuity  could  introduce  difficulty  in  the  application  of  FFT  techniques  because 
of  the  "wrap-around”  effect  inherent  to  the  FFT  process. 

If  an  end  discontinuity  is  present,  the  FFT  attempts  to  duplicate  the  time-domain 
discontinuity  by  introduction  of  excessively  large  high  frequency  components  to  the 
frequency  domain  characterization  of  the  wave  profile.  These  high  frequency  compo¬ 
nents  are  introduced  because  of  the  end  discontinuity,  but  are  propagated  through 
the  entire  wave  record.  After  application  of  the  Inverse  FFT,  the  resultant  time  trace 
has  the  spurious  high  frequency  content  spread  throughout  the  wave  history.  There 
is  no  damping  in  the  dispersion  relation  so  any  numerical  noise  introduced  by  end 
discontinuities  is  not  dissipated  as  it  is  propagated  to  other  spatial  locations  as  part 
of  the  wave  process. 

To  avoid  the  introduction  of  an  end  point  discontinuity,  a  portion  of  a  wave  cycle 
is  artificially  generated  by  the  program  in  this  Predict  option.  First,  the  period  of  a 
characteristic  wave  cycle,  Tchar,  is  calculated  by  averaging  the  wave  upcrossing  periods 
for  the  first  and  last  ten  cycles  in  the  wave  record.  How  large  a  portion  of  a  single 
characteristic  wave  cycle  is  required  for  continuity  is  determined  by  estimating  the 
“phases”  of  the  waves  at  the  beginning  and  ending  of  the  wave  record  as  ©i  and  ©2: 

©1  =  arctan  (^) 


and 


©2 


arctan 


(^) 
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Figure  4.2:  Phase  Calculation  for  Added  Portion  of  Wave  Cycle 


where 

X2  =  Average  offset 
Xi  =  Slope  between 
X2  —  Slope  between 


of  the  last  two  points  in  the  wave  record  from  the  mean 
the  first  two  point  in  the  wave  record 
the  last  two  points  in  the  wave  record 


The  portion  of  the  characteristic  wave  cycle  to  be  included  is  then  the  absolute 
value  of  (©1  —  ©a)  and  so  the  required  additional  time  is  given  by: 


(4.8) 


where  ^=21^  /Tchar 

The  required  number  of  zuiditional  time  points  is  then  Tadded!  (time  between  suc¬ 
cessive  time  points) 

A  third-order  polynomial  is  then  fit  through  the  last  two  points  in  the  wave  record 
such  that  the  newly  created  polynomial  would  be  continuous  with  the  first  two  points 
on  the  record  if  the  wave  record  repeated.  Thus,  a  small  number  of  artificial  data 
points  (necessarily  less  than  one  wave  cycle)  is  generated  which  makes  the  wave  record 
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continuous  between  its  end-points  so  no  spurious  high  frequency  energy  is  introduced 
into  the  FFT  and  IFFT  calculations. 

Note  that  when  predictions  are  made  for  locations  other  than  the  location  of  the 
original  wave  probe,  the  artificially  generated  wave  data  is  shifted  according  to  its 
Fourier  frequency  content,  as  is  the  rest  of  the  data.  Thus,  part  of  the  artificially 
generated  data  is  shifted  a  short  distance  into  the  body  of  the  actual  measured  data, 
from  which  it  generally  can  not  subsequently  be  removed. 


4.4.2  Simulate  Option 

End  discontinuities  are  not  problematic  in  the  Simulate  option  because  the  Simulate 
option  generates  waves  by  use  of  Fast  Fourier  Transform  techniques.  The  generated 
wave  components  have  the  convenient  property  that  the  resultant  time  trace  of  each 
of  the  component  waves  is  continuous  in  slope  and  offset  between  the  two  ends  of 
the  generated  time  series.  This  property  enables  the  generation  consistent  simulated 
waves  at  multiple  spatial  locations  without  requiring  special  treatment  of  end  discon¬ 
tinuities. 


4.4.3  Identify  Option 

End  discontinuities  could  be  problematic  in  the  Identify  option.  The  identification 
process  is  extremely  computationally  intensive.  The  input  wave  time  history  is  sep¬ 
arated  into  smaller  time  history  segments;  each  of  these  windows  generally  has  its 
length  set  to  an  even  power  of  two  to  reduce  the  total  computational  intensity  of 
the  FFT  process.  Discontinuities  will  generally  be  present  between  the  ends  of  these 
windows. 

A  different  approach  is  applied  to  gain  end-point  continuity  in  the  Identify  option 
such  that  the  length  of  each  window  can  be  maintained  as  an  even  power  of  two. 
The  approach  is  to  increase  the  length  of  each  window  by  introduction  of  artificially 
generated  points  on  each  end  of  the  each  wave  history  segment.  The  number  of 
artificial  data  points  plus  the  number  of  actual  wave  data  points  are  selected  such 
that  the  total  remains  an  even  power  of  two. 

The  artificial  data  is  generated  such  that  the  wave  frequency  of  the  end  of  the 
actual  wave  segment  is  maintained.  The  amplitude  of  the  artificial  waves  is  gradually 
reduced  to  zero,  so  the  window  has  the  appearance  of  a  segment  of  actual  wave  data 
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surrounded  at  each  end  by  a  ’’ramp”  of  artificial  data  which  reduces  each  end  of  the 
window  to  the  mean  of  that  window.  Continuity  between  the  ends  of  the  window 
is  assured  because  each  end  has  waves  of  zero  amplitude.  The  dispersion  relation  is 
not  applied  in  the  Identify  option,  so  dispersion  of  the  ramps  into  the  wave  record  is 
not  a  concern.  The  artificial  data  is  removed  from  the  ends  of  each  segment  prior  to 
concatenation  of  the  segments  into  the  identified  first  and  second  order  wave  histories. 
The  details  of  the  ramp  implementation  are  further  discussed  in  Chapter  3. 


4.5  Input  Specification 


The  input  specification  for  the  prediction  of  a  first-order  wave  process  is  in  a  command¬ 
line  format  similar  to  the  input  for  the  simulation  or  identification  options. 

Typical  input  file:  syntax  description 


#  Wave  Prediction  Input  File 

predict  filename  npts  dt 
define  gravity  value 
define  varlimit  value 
define  omgmax  value 
depth  value 

write  history  filenamel  filenameB 

location  n/oc 

valuel 

value2 

valuenloc 


#  Typical  input  file:  syntax  description 

Any  line  beginning  with  a  “#”  is  treated  as  a  comment  line  in  the  input  file  and  is 
ignored  by  the  program.  Blank  lines  are  also  ignored  by  the  program.  As  with  other 
options  of  Wavemaker,  most  lines  in  the  input  file  can  be  entered  in  any  sequence. 
There  are  two  exceptions:  the  first  non-comment  line  must  indicate  the  type  of  analy¬ 
sis  to  be  performed  (simulate,  identify  or  predict),  and  the  second  is  that  the  location 
specification  must  be  the  last  set  of  lines  in  the  file. 

predict  filename  npts  dt 

The  required  keyword  predict  indicates  the  user  intends  to  predict  the  wave  history  at 
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another  spatial  location  based  on  a  (measured)  wave  history.  This  command  requires 
three  sequential  arguments: 

•  filename:  A  character  string  indicating  the  name  of  the  file  containing 
the  observed  or  identified  wave  time  history  for  which  wave  predictions 
at  other  locations  are  to  be  performed.  The  data  in  the  first  column  in 
filename  is  the  wave  time  history.  Any  blank  lines  in  filename  or  lines  that 
do  not  begin  with  a  number  are  ignored. 

•  npts:  An  integer  number  indicating  the  number  of  data  lines  to  be  read 
from  the  input  wave  file,  npts  times  dt  equals  the  elapsed  time.  If  npts 
is  specified  to  be  a  number  smaller  than  the  number  of  lines  in  the  data 
file,  the  first  npts  data  points  in  the  file  are  used.. 

•  dt:  A  real  number  indicating  the  time  resolution  of  the  wave  history  pro¬ 
vided  in  filename,  dt  is  the  difference  between  measurement  times  of  two 
successive  data  points  in  the  input  wave  history. 


define 

The  optional  keyword  define  is  to  be  followed  by  another  keyword  indicating  which 
constant  is  to  be  defined. 

define  varlimit  value 
define  gravity  value 
define  omgmax  value 


The  keywords  varlimit,  gravity,  and  omgmax  have  the  same  meaning  as  in  the 
simulation  and  identification  sections.  The  value  of  omgmax  is  in  units  of  radians  per 
second  and  must  be  specified  identically  in  identify  and  in  predict  to  obtain  consistent 
results. 

depth  value 

The  required  keyword  depth,  as  in  the  simulation  option,  indicates  the  water  depth 
at  which  the  identification  analysis  is  to  be  performed,  value  is  a  real  number.  As 
in  the  simulate  and  identify  options,  depth  can  be  specified  in  any  units  consistent 
with  the  units  of  gravity  specified  using  the  define  keyword.  If  the  user  does  not 
explicitly  specify  a  value  for  gravity,  the  default  value  9.81  meters/sec^  and  then  the 
depth  value  should  be  specified  in  meters. 

write  history  filenamel  filename2 
write  statistics  filenames  filename4 
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The  optional  keyword  write  is  to  be  followed  by  either  the  keyword  history  or 
the  keyword  statistics,  write  is  used  to  specify  the  files  to  which  the  predicted  time 
histories  or  summary  statistics  are  to  be  written. 

write  history  filenamel  filename2 


The  final  length  of  each  of  these  files  should  be  the  same  as  the  length  of  the 
original  input  wave  history,  plus  a  small  amount  of  descriptive  text.  Output  file 
formats  are  covered  in  section  4.4  and  in  more  detail  in  2.4.1. 

filenamel  is  the  file  to  which  the  predicted  first-order  wave  histories  for  each  of  the 
specified  locations  are  written. 

filename2  is  the  file  to  which  the  predicted  second-order  wave  histories  for  each  of  the 
specified  locations  are  written. 

Default  values  are  assigned  to  filenamel  and  filename2  if  none  are  specified  by  the 
user.  These  default  filenames  are  gauss.his  and  ngauss.his,  respectively. 

write  statistics  filenames  filename^ 


The  final  length  of  each  of  these  files  should  be  the  same  as  the  number  of  ouput 
locations  specified  in  nloc  below,  plus  a  small  amount  of  descriptive  text.  Output  file 
formats  are  covered  in  section  4.4  and  in  more  detail  in  2.4.1. 

filenames  is  the  file  to  which  summary  statistics  of  the  predicted  first-order  wave 
histories  for  each  of  the  specified  locations  is  written. 

filename^,  is  the  file  to  which  summary  statistics  of  the  predicted  second-order  wave 
for  each  of  the  specified  locations  are  written. 

Default  values  are  assigned  to  filenames  and  filename4  if  none  are  specified  by  the 
user.  These  default  filenames  are  gauss.sta  and  ngauss.sta,  respectively. 

location  nloc: 

The  required  keyword  location  is  used  to  specify  the  number  of  spatial  locations  at 
which  both  the  first-  and  total  second-order  wave  histories  are  to  be  predicted. 

nloc  is  an  integer  specifying  the  desired  number  of  locations.  The  maximum 
allowable  number  of  locations  is  50. 
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value  1 
value2 

valuenloc 


Each  of  these  real  numbers  indicates  one  spatial  location  measured  in  the  along- 
wave  direction  at  which  prediction  results  are  to  be  generated.  Both  negative  and 
positive  values  are  acceptable.  Positive  values  are  down-stream  of  the  original  input 
wave  history;  negative  values  are  closer  to  the  wave  generating  device.  As  noted 
earlier,  this  specification  of  the  spatial  locations  must  be  after  all  other  input  in  the 
file. 


It  is  suggested  that  one  of  the  output  locations  be  0.0,  where  the  waves  were 
originally  measured.  It  can  then  be  confirmed  that  the  wave  data  has  been  correctly 
read  and  interpreted  by  the  program  by  verifying  that  the  predicted  wave  history  at 
the  wave  probe  equals  the  target  wave  history  measured  at  the  same  location. 


4.6  Output  Format 


A  total  of  four  output  files  are  produced  by  the  driver  program.  Two  output  files 
contain  the  time  histories:  one  for  the  underlying  first-order  wave  histories  and  the 
other  for  the  total  second-order  histories.  The  other  two  output  files  contain  wave 
statistics:  first  four  moments,  minimum  and  maximum.  Again,  results  for  the  first- 
and  second-order  wave  histories  are  separated  into  two  files. 


4.6.1  Time  History  Output 


As  noted  previously,  by  default  the  first-  and  second-order  histories  are  written  to  the 
files  gauss.his  and  ngauss.his.  Other  choices  of  output  filenames  can  be  specified 
by  the  optional  write  history  command.  The  format  of  this  output  depends  on  the 
number  of  spatial  locations  specified.  If  the  number  of  locations  (n/oc)  is  less  than  or 
equal  to  8  then  the  output  is  in  Formatl  otherwise  the  output  is  in  Format2.  Both 
of  these  formats  write  out  3  header  lines  beginning  with  a  sign.  These  are  to  be 
treated  as  comment  lines  in  the  output  file. 

Formatl  outputs  data  in  n/oc-t-1  columns.  The  length  of  each  column  is  equal 
to  the  number  of  points  desired  in  each  prediction.  The  first  column  contains  the 


44  Chapter  4.  Prediction  of  Second-Order  Random  Waves 

time  increments  in  seconds  going  from  0  to  T  with  dt  =  T/npts.  Columns  2  through 
n/oc+1  contain  the  predicted  time  history  values  at  the  specified  locations  xl,  x2,  x3, 
... xnloc,  respectively.  Thus,  column  2  contains  the  wave  elevation  at  location  xl, 
column  3  contains  wave  elevation  at  location  x2,  and  so  on. 

Format2  is  for  handling  nloc  greater  than  8.  The  output  begins  with  the  time 
increment  Ti  in  seconds  on  a  line  by  itself.  The  time  history  values  for  the  specified 
spatial  locations  at  time  Ti  are  written  in  the  next  line  onwards,  in  sets  of  10.  So  if  9 
locations  were  specified  (i.e.,  nloc  =  9)  then  the  time  increment  is  printed  on  a  line  by 
itself  followed  by  a  line  containing  9  time  history  values  at  that  time  increment.  The 
next  line  contains  the  next  time  increment  followed  by  another  set  of  9  values,  and 
so  on.  If,  on  the  other  hand  say  28  locations  were  specified,  then  a  time  increment  is 
written  on  a  line  followed  by  28  time  history  values  (corresponding  to  28  locations  at 
that  time  increment)  in  the  next  3  lines.  The  first  line  of  the  3  lines  contains  10  time 
history  values  for  the  first  10  locations  specified.  The  next  line  contains  10  history 
values  for  locations  11  through  20  and  the  following  line  which  is  the  third  line  of  the 
set  will  contain  only  8  history  values  for  location  21  through  28. 


4.6.2  Wave  Statistics  Output 

The  statistics  of  the  simulated  histories  are  also  estimated  by  the  driver  program. 
These  statistics  include  the  mean,  standard  deviation,  skewness,  kurtosis,  minimum, 
and  maximum.  As  noted  in  the  previous  section,  first-  and  second-order  simulation 
results  are  written  by  default  to  the  files  gauss.sta  and  ngauss.sta,  respectively. 
The  optional  command  write  statistics  can  alter  this  choice  of  output  filenames. 

The  output  format  in  both  of  these  files  begins  with  2  header  lines,  each  of  which 
begins  with  a  sign.  The  output  is  in  seven  columns.  The  first  column  specifies 
the  spatial  location.  The  following  six  columns  contain  statistics  of  the  wave  history 
at  the  spatial  location  specified  in  column  1.  Columns  2  through  7  contain,  the  mean, 
standard  deviation,  skewness,  kurtosis,  minimum, ,  and  maximum  value  in  that  order. 


4.7  Examples 


In  this  section  we  present  two  sample  problems  to  illustrate  the  use  of  the  prediction 
capabilities  of  wavemaker.  Example  1  is  based  on  the  examples  presented  in  the  sim¬ 
ulation  and  identification  chapters.  Sample  input  and  output  files  of  this  prediction 
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example  are  included  on  the  distribution  diskette.  Example  2  presented  here  demon¬ 
strates  the  prediction  of  first-  and  second-order  components  of  a  measured  wave  tank 
history.  Note  that  sample  input  or  output  files  of  this  second  example  are  not  included 
in  the  distribution. 


4.7.1  Example  1 

The  example  of  the  simulation  and  identification  capabilities  of  wavemaker  involved 
simulating  a  second-  order  wave  history  characterized  by  a  JONSWAP  spectrum  with 
H,  =  12  [m],  Tp  =  14s  and  7  =  3.3  in  a  70  [m]  water  depth.  We  will  use  the  first-order 
wave  history  and  try  to  predict  the  first-  and  second-order  wave  history  at  another 
location  for  comparison  with  the  histories  directly  simulated  at  the  same  alternate 
location  (60  meters  down-stream). 

The  input  file  for  the  prediction  run  is: 


#  Wave  Prediction  Input  File 

predict  gauss.ide  4096  0.5 
depth  70.0 

write  history  gauss.pre  ngauss.pre 
define  omgmax  3.0 
location  3 
0.0 
2.0 
60.0 


The  input  file  gauss.ide  contains  a  column  of  4096  real  numbers  (see  sample  files 
listed  in  the  appendix)  having  an  incremental  time  step  of  0.5  seconds.  This  file  is 
read  in  as  the  first  order  components  of  the  identified  observed  wave  history. 

The  linear  (gaussian)  prediction  are  written  to  the  file  gauss.his,  while  the  predic¬ 
tion  including  the  non-linear  (non-gaussian)  terms  is  written  to  the  file  ngauss.his. 

Figure  3.2  shows  a  portion  of  the  observed  wave  history  spectrum  and  correspond¬ 
ing  first-  and  second-order  wave  histories.  We  see  that  small  second-order  contribution 
pointed  out  in  the  identify  example  contributes  meaningfully  to  the  maximum  crest 
heights. 
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Time  (sec) 

Figure  4.3:  Non-Linear  Prediction  and  Simulation  Results 


4.7.2  Example  2 


In  this  example  we  will  continue  example  2  from  the  Identify  section  by  predicting 
the  wave  profile  at  two  locations  in  addition  to  the  wave  measurement  location.  The 
underlying  first-order  wave  component  for  a  measured  wave  tank  history.  For  this 
example  the  identified  history  is  located  in  file  gauss.ide  and  has  a  dt  =  0.3354 
seconds.  We  will  use  windows  of  winsize  =  512  to  identify  the  first-order  components. 
The  input  to  wavemaker  is: 


Elevation  (meters) 
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Figure  4.4:  Second-Order  Predicted  Results  -  Wave  Tank  Example 


#  Wave  prediction  Input  File 

predict  gauss.ide  4096  0.3354 
depth  175.0 
define  omgmax  4.2 
write  history  gauss2.pre  ngauss2.pre 
location  3 
0.0 
2.0 
60.0 


Figure  4.4  shows  a  portion  of  the  wave  history  as  predicted  at  three  special  lo¬ 
cations.  Note  that  the  shape  of  the  wave  profile  has  changed  only  slightly  when  it 
is  shifted  a  short  distance  (2  meters),  but  has  more  significant  shape  changes  at  the 
larger  shift  (60  meters). 
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Distribution 


The  WAVEMAKER  routine  and  example  files  have  been  distributed  on  a  DOS  formatted 
3.5  inch  floppy  diskette.  The  diskette  contains  the  source  code  files  (of  the  form  ♦.f), 
example  input  (♦.inp),  and  output  files  (*.sta,  *.his,  and  ♦.ide). 


5.1  Copying  the  Diskette 


Copy  the  contents  of  the  diskette  on  to  your  host  computer  (computer  on  which  you 
will  run  WAVEMAKER).  After  the  copying  is  done,  your  host  computer  should  have: 

•  Example  input  files  identify .  inp,  predict .  inp,  simulate .  inp  and  hist  .dat. 

•  Example  output  files  gauss. his,  ngauss.his,  gauss. sta,  ngauss.sta, 
gauss.ide,  ngauss.ide,  ngauss.pre,  gauss. pre,  identify. out,  predict  .out, 
and  simulate . out. 

•  all  the  source  files  *.f  and  Makefile 

Table  5.1  shows  the  input  and  output  files  that  are  specific  to  the  simulation,  identi¬ 
fication  and  prediction  examples. 


5.2  Compiling  the  Source 


On  a  Unix  workstation,  the  Makefile  can  be  used  to  build  the  WAVEMAKER  executable. 
To  compile  on  your  host  Computer; 
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Table  5.1:  Distributed  Files  for  Simulation  and  Identification  Examples 


Simulation  Example 

File  Type 

File  Name 

Description 

Input  File 

simulate . inp 

Input  to  WAVEMAKER 

Output  Files 

gauss. his 

Simulated  first-order  wave  histo¬ 
ries  at  each  specified  location 

ngauss.his 

Simulated  combined  first-  and 
second-order  wave  histories  at 
each  specified  location 

gauss. St a 

Statistics  of  the  simulated  first- 
order  histories  at  each  specified  lo¬ 
cation 

ngauss . sta 

Statistics  of  the  simulated  com¬ 
bined  histories  at  each  specified  lo¬ 
cation 

Identification  Example 

File  Type 

File  Name 

Description 

Input  Files 

identify . inp 

Input  to  WAVEMAKER 

hist.dat 

Observed  wave  history  for  which 
underlying  first-order  history  is  to 
be  identified 

Output  Files 

gauss . ide 

Identified  first-order  wave  history 

ngauss . ide 

Identified  second-order,  combined 
first-  and  second-order,  and  ob¬ 
served  wave  histories 

Prediction  Example 

File  Type 

File  Name 

Description 

Input  Files 

predict . inp 

Input  to  WAVEMAKER 

Output  Files 

gauss. pre 

Predicted  first-order  wave  history 

ngauss. pre 

Predicted  combined  first-  and 
second-order  wave  histories  at 
each  specified  location 

5.3.  Executing  the  Routine 
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•  change  directory  to  the  subdirectory  containing  the  source  code  files 

•  type  “make  wavmkr”  (without  the  double  quotes)  at  the  Unix  prompt 
and  press  return 

The  above  will  compile  all  the  files  listed  in  Makefile  and  link  them  to  make  the 
executable  wavmkr.  Note  that  the  executable  wavmkr  is  still  in  the  current  directory, 
and  may  be  moved  to  the  directory  which  contains  the  input  file,  or  you  can  specify 
the  path  to  the  executable  file  in  order  to  run  it  from  other  directories. 

On  other  operating  systems  and  architectures,  follow  whatever  is  the  standard 
procedure  for  compiling  and  linking  FORTRAN  source  code  distributed  over  multiple 
files. 

WAVEMAKER  has  been  developed  on  a  Sun  Sparcstation2,  using  version  1.4  of  the 
Sun  FORTRAN  compiler.  Every  attempt  has  been  made  to  adhere  to  the  ANSI 
FORTRAN  77  standard,  ensuring  portability  of  the  code. 


5.3  Executing  the  Routine 


At  the  Unix  prompt  type  “wavmkr  <  simulate.inp”  (without  the  double  quotes) 
in  order  to  execute  Wavemker  and  perform  the  simulation  example  analysis.  The 
program  reads  input  from  the  standard  logical  input  unit.  The  logical  units  for 
standard  input,  standard  output,  and  standard  error  are  all  used  in  WAVEMAKER.  On 
the  Sun  compiler,  0  is  used  for  standard  error,  5  is  standard  input  and  6  is  standard 
output.  If  the  appropriate  unit  numbers  are  different,  they  can  be  set  using  the  lOER, 
lOIN,  and  lOOU  variables  in  the  WAVEMAKER  driver  program,  and  the  package  can  be 
recompiled. 

Run  the  example  problem  using  the  complied  code  to  check  if  you  get  the  same 
output  as  provided  in  the  example  output  files.  Note  that  wavmkr  will  overwrite  any 
existing  file  if  you  specify  its  name  as  the  target  output  file  using  the  write  command 
in  the  input  file. 

Similarly,  at  the  Unix  prompt  type  “wavmkr  <  identify.inp”  (without  the  double 
quotes)  in  order  to  execute  WAVEMAKER  and  perform  the  identification  example  anal¬ 
ysis.  The  resulting  output  histories  can  be  compared  to  the  corresponding  example 
output  histories  to  verify  successful  compilation  of  WAVEMAKER. 

Finally,  at  the  Unix  prompt  type  ‘Svavmkr  <  predict.inp”  (without  the  double 
quotes)  in  order  to  execute  WAVEMAKER  and  perform  the  prediction  example  analysis. 
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The  resulting  output  histories  can  again  be  compared  to  the  corresponding  example 
output  histories  to  verify  successful  compilation  of  WAVEMAKER. 
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Appendix  A.  Output  Files  for  Simulation  Example 


Output  File:  gauss. his  _ 

#  Underlying  First-Order  Wave  Process 

#  Wave  Elevation  at  Spatial  Location 


Time (sec. ) 

0.00 

60.00 

0.000000 

-3.189 

-4.754 

0.500000 

-1.754 

-5.466 

1.000000 

-0.239 

-5.410 

1.500000 

0.906 

-4.667 

2.000000 

1.436 

-3.660 

2.500000 

1.681 

-2.714 

3.000000 

1.623 

-1.676 

3.500000 

1.595 

-0.875 

4.000000 

1.862 

-0.137 

4.500000 

2.282 

0.582 

5.000000 

2.802 

1.211 

5.500000 

3.230 

1.379 

e.  000000 

3.263 

1.263 

6.500000 

2.448 

1.252 

7.000000 

1.403 

1.572 

7.500000 

0.553 

2.166 

8.000000 

-0.121 

2.447 

8.500000 

-0.773 

3.076 

9.000000 

-1.216 

3.472 

9.500000 

-1.588 

3.272 

10.000000 

-2.433 

2.507 

10.500000 

-3.543 

1.492 

11.000000 

-4.310 

0.022 

11.500000 

-4.343 

-1.441 

2046.500000 

-3.871 

0.056 

2047.000000 

-4.266 

-1.729 

2047.500000 

-4.040 

-3.447 

Note  that  indicates  more  numbers.  Since  the  files  are  long,  we  present  truncated  versions 
of  the  output  files. 
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Output  File:  ngauss.his 


#  Total  Second- 

-Order  Wave  Process 

# 

Wave  Elevation  at  Spatial  Location 

#  Time (sec.) 

0.00 

60.00 

0.000000 

-3.561 

-4.774 

0.500000 

-2.670 

-4.904 

1.000000 

-0,945 

-4.906 

1.500000 

0.892 

-4.458 

2.000000 

1.894 

-3.636 

2.500000 

2.154 

-2.719 

3.000000 

1.949 

-1.632 

3.500000 

1.835 

-0.810 

4.000000 

2.029 

-0.150 

4.500000 

2.394 

0.494 

5.000000 

2.982 

1.263 

5.500000 

3.504 

1.740 

6.000000 

3.399 

1.748 

6.500000 

2.247 

1.508 

7.000000 

1.054 

1.561 

7.500000 

0.447 

2.191 

8.000000 

0.119 

2.688 

8.500000 

-0.431 

3.386 

9.000000 

-0.983 

3.711 

9.500000 

-1.462 

3.471 

10.000000 

-2.379 

2.563 

10.500000 

-3.558 

1.095 

11.000000 

-4,273 

-0.814 

11.500000 

-4.100 

-2.098 

12.000000 

-3.351 

-2.355 

12.500000 

-2.713 

-2.254 

13.000000 

-2.058 

-2.384 

13.500000 

-0.726 

-2.484 

2046.000000 

-3.361 

1.599 

2046.500000 

-3.657 

-0.860 

2047.000000 

-3.765 

-3.229 

2047.500000 

-3.710 

-4.462 
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Appendix  B.  Input/Output  Files  for  Identification  Example 


Note  that  indicates  more  numbers.  Since  the  files  are  long,  we  present  truncated  versions 
of  the  output  files. 


Output  File:  gauss.ide _ 

#  Indentified  First-Order  Wave  History 
-2.5860 
-1.2615 
-0.0398 
0.7557 
1.0355 
1.2597 
1.3753 
1.5219 
1.8287 
2,2055 
2.7083 
3.1522 
3.1924 
2.3991 
1.4205 
0.6076 
-0.1573 
-0.9790 
-1.5084 


-0.9156 
-1.4972 
-1.6011 
-1.5946 
-1.9273 
-2.3777 
-2; 9893 
-3.6520 
-4.3427 
-4.2629 
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Output  File:  ngauss.ide 


#  Identified  2nd~0rder 

Wave,  Identified  Total  Wave,  Input  Wave  History 

-0.6029 

-3.1890 

-3.1890 

-0.4924 

-1.7539 

-1.7540 

-0.1992 

-0.2390 

-0.2390 

0.1502 

0.9060 

0.9060 

0.4005 

1.4359 

1.4360 

0.4212 

1.6810 

1.6810 

0.2477 

1.6230 

1.6230 

0.0731 

1.5950 

1.5950 

0.0333 

1.8620 

1.8620 

0.0765 

2.2820 

2.2820 

0.0937 

2.8020 

2.8020 

0.0778 

3.2300 

3.2300 

0.0706 

3.2630 

3.2630 

0.0489 

2.4480 

2.4480 

-0.0175 

1.4030 

1.4030 

-0.0546 

0.5530 

0.5530 

0.0362 

-0.1210 

-0.1210 

0.2060 

-0.7730 

-0.7730 

0.2924 

-1.2160 

-1.2160 

-0.1915 

-1.1071 

-1.1070 

0.3880 

-1.1092 

-1.1090 

0.5239 

-1.0772 

-1.0770 

0.3025 

-1.2921 

-1.2920 

-0.0245 

-1.9518 

-1.9520 

-0.2870 

-2.6647 

-2.6650 

-0.3756 

-3.3649 

-3.3650 

-0.2191 

-3.8711 

-3.8710 

0.0766 

-4.2661 

-4.2660 

0.2229 

-4.0400 

-4.0400 
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Appendix  C.  Output  Files  for  Prediction  Example 


Output  File:  gauss.pre 


-X"  -  * 

#  Underlying  First-Order 

Wave  Process 

# 

Wave  Elevation  at 

Spatial  Location 

#  Time (sec.) 

0.00 

2.00 

60.00 

0.000000 

-2.586 

-2.764 

-4.224 

0.500000 

-1.262 

-1.813 

-3.682 

1.000000 

-0.040 

-0.522 

-3.141 

1.500000 

0.756 

0.631 

-2.784 

2.000000 

1.036 

1.271 

-2.384 

2.500000 

1.260 

1.380 

-1.878 

3.000000 

1.375 

1.431 

-1.262 

3.500000 

1.522 

1.560 

-0.877 

4.000000 

1.829 

1.782 

-0.709 

4.500000 

2.205 

2.023 

-0.575 

5.000000 

2.708 

2.499 

-0.244 

5.500000 

3.152 

2.819 

0.405 

6.000000 

3.192 

3.067 

0.899 

6.500000 

2.399 

2.804 

1.887 

7.000000 

1.421 

1.963 

2.813 

7.500000 

0.608 

0.738 

3.390 

8.000000 

-0.157 

-0.130 

3.582 

2043.000000 

-0.916 

-0.800 

5.185 

2043.500000 

-1.497 

-1.539 

5.590 

2044.000000 

-1.601 

-1.830 

5.317 

2044.500000 

-1.595 

-1.879 

4.422 

2045.000000 

-1.927 

-1.921 

3.287 

2045.500000 

-2.378 

-2.146 

2.461 

2046.000000 

-2.989 

-2.848 

1.366 

2046.500000 

-3.652 

-3.387 

-0.389 

2047.000000 

-4.343 

-4.012 

-2.238 

2047.500000 

-4.263 

-4.401 

-3.703 
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Output  File:  ngauss.pre 


#  Total  Second-Order  Process 

#  Wave  Elevation  at  Spatial  Location 


#  Time (sec.) 

0.00 

2.00 

60.00 

0.000000 

-3.205 

-2.977 

-3.840 

0.500000 

-2.255 

-2.572 

-3.420 

1.000000 

-0.736 

-1.479 

-2.968 

1.500000 

0.777 

0.088 

-2.594 

2.000000 

1.522 

1.442 

-2.279 

2.500000 

1.667 

1.953 

-2.001 

3.000000 

1.498 

1.863 

-1.437 

3.500000 

1.550 

1.652 

-0.786 

4.000000 

1.914 

1.738 

-0.370 

4.500000 

2.278 

2.054 

-0.316 

5.000000 

2.713 

2.616 

-0.246 

5.500000 

3.195 

2.978 

0.220 

6.000000 

3.338 

3.253 

0.576 

6.500000 

2.522 

2.923 

1.419 

7.000000 

1.401 

1.882 

2.427 

7.500000 

0.544 

0.539 

3.426 

8.000000 

-0.073 

-0.151 

3.971 

2043.000000  -1.124  -1.700  5.807 


2043.500000 

-1.049 

-1.599 

5.835 

2044.000000 

-0.994 

-1.140 

4.913 

2044.500000 

-1.239 

-1.159 

3.998 

2045.000000 

-1 . 922 

-1.718 

3.418 

2045.500000 

-2.607 

-2.372 

2.621 

2046.000000 

-3.261 

-3.101 

0.510 

2046.500000 

-3.781 

-3.454 

-2.323 

2047.000000 

-4.255 

-3.946 

-4.103 

2047.500000 

-4.025 

-4.288 

-4.506 
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